A correlator product state study of molecular magnetism 
in the giant Keplerate Mo72Fe3o 
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We have studied the properties of the giant Keplerate molecular magnet Mo72Fe3o, as a function of 
applied magnetic field, using the correlator product state (CPS) tensor network ansatz. The magnet 
is modeled with an S = 5/2 antiferromagnetic Heisenberg Hamiltonian on the 30-site icosidodec- 
ahedron lattice, a model for which exact diagonalization is infeasible. The CPS ansatz produces 
significant improvements in variational energies relative to previous studies using the density matrix 
renormalization group, a result of its superior ability to handle strong correlation in two dimensional 
spin systems. The CPS results reaffirm that the ground state energies adhere qualitatively to the 
parabolic progression of the rotational band model (RBM), but show important deviations near 1/3 
of the saturation field. These deviations predict anomalous behavior in the differential magnetiza- 
tion and heat capacity that cannot be explained by the RBM alone. Finally, we show that these 
energetic deviations originate from a qualitative change in the ground state that resembles a finite 
size analogue of a phase transition. 



I. INTRODUCTION 



Molecular magnets are classic examples of chemical 
systems containing a large number of localized, strongly 
correlated electrons. Their study has been motivated 
both from potential applications in storage and quantum 
computing, as well as by the fundamental challenges as- 
sociated with their chemical synthesis and their physical 
magnetic properties [THS]. In recent years, using polyox- 
omctalate chemistry [7] , some very large molecular mag- 
nets have been synthesized [51 H] . These so-called giant 
Keplerate magnets earn their name from the geometric 
arrangement of the ions, which lie at the vertices of reg- 
ular solids. The largest such magnet made to date is 
based on the icosidodecahedron, and consists of corner 
sharing triangles arranged around pentagons (see Figure 
fl). The metal species can be varied, and magnets in- 
cluding V, Cr, and Fe ions have been made, although the 
Fe based Keplerate magnet has been the most studied so 
far pHUlS) . The corner sharing triangle geometry leads 
to magnetic frustration and unusual magnetic properties 
[B] which are of interest in this work. 

The theoretical description of magnetism in the Ke- 
plerate magnets is extremely challenging. The basic rea- 
son is the size of the Hilbert space associated with the 
magnetic centers. In the case of the Feao-Keplerate, each 
Fe center is a 3-1- ion with 5 unpaired spins in a near 
perfect octahedral coordination, and we can view each 
center as effectively e.n S — 5/2 spin [12]. Arranging 
the spins on the vertices of the icosidodecahedron (see 
Figures [I] and ^ , we model their interactions using the 
Heisenberg Hamiltonian, 



H — J y Si ■ Sj 



(1) 




where {ij) represents a summation over nearest neigh- 
bors. Since there are 30 S' = 5/2 spins, the correspond- 
ing Hilbert space is of dimension 6'^", or roughly 10^"^, a 
mole of quantum states! This is far too large to treat 



FIG. 1: The giant Keplerate Mo72Fe3o molecular magnet is 
shaped lilce an icosidodecahedron, with the Fe atoms posi- 
tioned on the vertices and the — O — Mo — O— bridges along 
the edges. Three example correlators have been shaded in 
red: a two-site nearest neighbor, a three-site triangle, and a 
five-site bow tie. 



using the exact diagonalization methods that are usually 
employed for molecular magnets [19j . 

In this work, we use a variational methodology based 
on correlator product states (CPS) |20U21j . in conjunction 
with the Heisenberg Hamiltonian, to model the low-lying 
states of the Fcao-Keplerate magnet. Correlator product 
states, also known as entangled plaquette states :22;, 23] 
or complete graph tensor networks [23] , provide a simple 
approximation to the quantum wavefunction amplitude 
for a large number of spins as a product of amplitudes of 
snraller overlapping subsets of spins. The term correla- 
tor refers to the amplitudes on the subsets of spins. The 



CPS approximation derives from an attempt to gener- 
alize the density matrix renormalization group (DMRG) 
pS] , a powerful method for strongly correlated electrons 
that has been applied both to realistic quantum chemical 
problems [2SMT] as well as many model condensed mat- 
ter Hamiltonians [32H34J. Exler and Schnack previously 
used the DMRG to study the Feso-Keplerate magnet [13] , 
providing a qualitative demonstration of the existence of 
a quantum rotational band. However, the DMRG has 
difhculty in accurately treating large systems where cor- 
relations are not ordered in a one-dimensional fashion. 
Unlike the DMRG, the CPS is not biased towards one- 
dimensional correlations, and thus in principle can be an 
efficient ansatz for the correlations present in the Kepler- 
ate magnets. Here we will compare our CPS calculations 
not only to the available experimental measurements, but 
also to the earlier theoretical DMRG work of Exler and 
Schnack, demonstrating the improved ability of the CPS 
to describe correlations in general systems. 

The structure of our study is as follows. We first give 
an overview of the theoretical and experimental results 
for magnetism in the giant Keplerate magnets and dis- 
cuss, in particular, features related to magnetic frustra- 
tion in the icosidodecahedron (Section [h]). We then de- 
scribe the general theory behind the CPS wavefunction 
(Section IH A I , how it is optimized via variational Monte 



tion behavior (Section IV D). Finally, we conclude with 



Carlo (Section IIIBl, and the specific form of the wave- 



function we use in this work (Section HIC). We next 



present the results of our calculations in light of exper- 
imental and earlier theoretical work on the magnet. In 
particular, we present total energies (Section IV A I, low 
temperature properties (Section IV B ) , spin correlations 
(Section IV C ) , and an analysis of possible phase transi- 




FIG. 2: The icosidodecahedron lattice flattened to a planar 
graph. The vertices are shown under the particular three- 
coloring that was used to generate the initial guess for our 
CPS wavefunction optimization. 



some perspectives for further work on the Keplerate sys- 
tems, and ways to generalize the CPS approach to other 
complex molecular systems (Section [v|. 



II. MAGNETISM IN THE GIANT 
KEPLERATES 

Keplerate systems are interesting from the viewpoint 
of quantum magnetism due to the presence of frustration 
effects [B]. One way to define a frustrated magnet is 
one where the classical Ising model, whose spins only 
assume up (u) and down (d) orientations, has a large 
degeneracy. This is the case for triangles, where the uud, 
udu, and duu configurations are all degenerate. In the 
classical Heisenberg model, where spins can point in any 
orientation, the spin triangle has a continuous manifold of 
degenerate ground states. In these states, the three spins 
are coplanar and rotated 120° from each other, and it is 
the orientation of the plane that creates the continuous 
degeneracy. 

While a single spin triangle already shows some frus- 
tration effects, such effects become even more pro- 
nounced in the case of corner sharing triangles [S] . This 
is the motif underlying the icosidodecahedron, whose sur- 
face consists of corner sharing triangles arranged around 
pentagons. In fact, the icosidodecahedron is the largest 
member of a family of Platonic solids, which also in- 
cludes the cuboctahedron and the truncated tetrahedron, 
whose surfaces are built from corner sharing triangles. 
The quantum Ising model on these lattices is highly frus- 
trated. These zero dimensional systems are especially 
important as they exist as finite size surrogates for their 
bulk planar counterparts, such as the two dimensional 
Kagome lattice [35 , which are believed to underlie ex- 
otic magnetism in solids. Because of their small size, the 
Platonic solid models allow the effects of corner sharing 
triangle frustration to be studied in an experimentally re- 
alizable system that is also accessible to many theoretical 
approaches. 

We now give a brief overview of some of the interesting 
properties that can arise from spin frustration in corner 
sharing triangle systems. One class of frustration effects 
is the presence of anomalies that occur at applied mag- 
netic field strengths close to 1/3 of the saturation field 
Bsat- (The saturation field is the field strength above 
which the ground-state has all spins aligned with the 
field) . It has been observed both experimentally [T71 |TS] 
and theoretically [171 Ull ES] that the differential suscep- 
tibility dM/dB (the rate of change of the total system 
magnetization with respect to field strength) displays a 
depression near field strengths of Bgat/i- A rough un- 
derstanding of this is that near Bsat/S^ the magnetically 
stiff uud states of the spin triangles become energetically 
competitive with the usual ground-state, but this alone 
only gives a qualitative accounting of the experimental 
data. In Keplerate systems, a more quantitative match 
T8" to the observed dM/dB depression was achieved un- 
der the assumption of random variations in the spin cou- 



plings within the classical Heisenberg model. Note that 
the Bsat/S anomaly does not only appear in the differen- 
tial susceptibility, but also shows up, for example, in zero 
temperature magnetization predictions |12| . in the heat 
capacity [3B], and as a phase transition in the classical 
Heisenberg model on the Kagome lattice [36l |37] . 

Another interesting aspect of frustrated spin systems 
is the possibility of unusually low-lying singlet excited 
states. Although these states will not be treated in this 
study, they have attracted a great deal of interest [121 
1161 I38j and are implicated as a means to explain puzzling 
experimental neutron scattering data jlOl 111) in the giant 
Keplerate magnets. 

To a first approximation (although see Ref. [15] ) mag- 
netism in these systems can be described by an isotropic 
Heisenberg model, with the M ion coupled antiferromag- 
netically via the Mo-0 bridges. The M ions are believed 
to lie in near perfect 0^ coordination with the oxygens, 
and the V, Cr, and Fe giant Keplerates can be thought of 
as S* = 1/2, 3/2, and 5/2 spin centers. Most experimen- 
tal work has focused on characterizing the S — 5/2 Fe 
system, and it is the corresponding S = 5/2 Heisenberg 
model on the icosidodecahedron to which we apply the 
CPS wavefunction. 

The essential problem in studying the icosidodecahe- 
dron Heisenberg model is the very large Hilbert space 
that needs to be considered, which is 2'^'^, 4^°, 6^° for 
the V, Cr, Fe species. Exact (full) diagonahzation of 
the Heisenberg Hamiltonian has been carried out in the 
5=1/2 case of the V30 magnet but is impossible for 
the other magnets. Nonetheless, many of the qualitative 
features of these systems appear to be well described by 
a rather simple model known as the quantum rotational 
band model. The Keplerate magnets are tripartite (see 
Figure [2]), and we can consider therefore a family of spins 
living on the A, B, C sub-lattices. The quantum rota- 
tional band model (REM) ,14j asserts that the energies 
of the states can be modeled as arising from the couplings 
of total spins on the A and B and C lattice as an effective 
triangle and is given by the Hamiltonian 



H, 



band 



= J- 



N 



s'~j(sl + sl + sl 



(2) 



where N is the number of spins, D and 7 are free param- 
eters, S is the net lattice spin, and Sa, Sb, and Sc are 
the net spins on each sublattice. The eigenstates of this 
Hamiltonian have energies 



£'(5', Sa, Sb,Sc) 
D 



J- 



N 



5(5 + 1) -7 I Y. Sq{Sq + 1) 

Qe{A.,B,C} 



(3) 



with degeneracies given by the number of ways a given 
total spin 5 can be made up from the sublattice spins 
via spin coupling rules. As will be discussed in the re- 
sults below, the REM has been successful at qualitatively 
reproducing experimental magnetizations |14j as well as 
the total energies of some higher level theoretical treat- 
ments ITS]. However, as we will demonstrate in the case 



of differential susceptibility, it fails to predict the pecu- 
liar properties of the Mo72Fe3o Keplerate magnet related 
to the i?sat/3 anomaly. For these effects we need to con- 
sider all the spin degrees of freedom, for which we need 
explicit approximations for the quantum wavefunction. 



III. CORRELATOR PRODUCT STATES 

A. General Theory 

Consider a set of k spins si . . . Sk- In an 5 = 5/2 sys- 
tem, such as in the Feso-Keplerate magnet, each s varies 
over the 6 mg levels of each iron center. The quantum 
wavefunction written in the complete spin Hilbert space 

is 

l*> = Yl *«1^2...«JS1S2-Sfc) (4) 

SiS2...Sfc 

s 

where s denotes the vector of spin configurations 

S1S2 ...Sk- 

The amplitude ^siS2...Sfc is infeasible to obtain exactly 
for a system as large as the Feao-Keplerate magnet. Cor- 
relator product states provide an approximation for the 
full amplitude in terms of simpler objects known as cor- 
relators. In spin systems, a correlator defines a set of 
amplitudes over a subset (domain) of the spin sites. For 
example, a correlator on sites i, j defines a set of ampli- 
tudes Cg-g.. Correlators can be constructed to act on an 
arbitrary number of sites (see Figure [I]) . Such a general 
correlator is written as Cs^ where s^ denotes the spin con- 
figuration of the subset of sites A. To obtain the CPS, we 
approximate the wavefunction amplitudes 'I's in Eq. Q 
as a product of correlator amplitudes over the different 
subsets of sites A, 



*« 



H' 



(5) 



Note that the domains A of the different correlators will 
usually contain overlapping sites. For example, a CPS 
wavefunction for a one-dimensional arrangement of spins 
with "nearest neighbor" correlators, would be written as 



*,, 



(6) 



Ey using correlators that cover increasingly larger num- 
bers of sites, we can make the CPS approximation arbi- 
trarily exact. 



B. Monte Carlo Optimization 

We use the variational Monte Carlo algorithm to op- 
timize the CPS wavefunction to obtain approximate 
ground-states of the Keplerate magnet. (We have shown 
elsewhere that the CPS wavefunction can also be used 
with non-stochastic algorithms [3T|, although these are 



not employed here). In variational Monte Carlo, the en- 
ergy is written as 



E 






(7) 



E 



(*i*; 



^l(s) 



where the local energy £^l(s) is defined by 



*s 



^l(«)-E^(«i^n')- 



(8) 



As long as ^s can be evaluated efficiently, which is the 
case for the CPS wavef unctions, a Markov chain can be 
used to sample the probability distribution |^sP/(^|^) 
and efficiently compute the overall energy as an aver- 
age of the sampled local energies. The energy is then 
variationally minimized using stochastic estimates for the 
gradient with respect to the correlator amplitudes. Note 
that it is easy to constrain the Monte Carlo sampling 
over s, for example in Eq. ([t]), to only those configu- 
rations with a given value of 5'^, and this allows us to 
obtain approximate ground-states in different Sz sectors. 
Once the wavefunctions are obtained, expectation val- 
ues for various correlation functions can also be readily 
computed by Monte Carlo sampling. 
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FIG. 3: Ground state energies for the S = 5/2 Heisenberg 
model on the icosidodecahedron for different total Sz sectors 
of the Hilbert space. In the main panel the CPS wavefunc- 
tion's energies are shown along with the corresponding fit to 
the RBM form, Eq. (|3|. In the inset the CPS energies are 
compared to the RBM produced by fitting DMRG energies 
[in]. See Section [TVAl 



IV. RESULTS 



C. Wavefunction and Optimization Details 



A. Total energies 



To study the FesQ-Keplerate magnet we used a CPS 
in the form of Eq. ((sl) with bow tie shaped correlators. 
There are 30 different bow ties in all, each defined by 
choosing one site and all of its nearest neighbors (see 
Figure [l] for an example) . After discovering that ran- 
domly chosen correlator amplitudes were not effective 
as initial guesses for the variational optimization, we 
chose to use as our guess a relatively simple state simi- 
lar to the classical ground state. To be precise, our ini- 
tial guess for Sz — ^ was chosen to be a spin-coherent 
state [39] which can be exactly represented by a CPS. In 
a spin-coherent state, the wavefunction amplitude fac- 
torizes into a product of amplitudes on individual sites 
^siS2.-Sfc = ^sx^^s'i • • ■ Csfc, and each site amplitude Cg^ de- 
fines a direction for the spin on the site. Here we chose 
the rotation angles for each site to be the classical ground 
state's spin direction for that site's sublattice (the sub- 
lattices were assigned based on the coloring shown in 
Figure pi). Starting from this guess, we optimized the 
wavefunction's energy under 5*2 = projection, and then 
used the resulting wavefunction as an initial guess for the 
Sz = ^ sector. In this fashion we worked our way up the 
magnetization ladder, obtaining a wavefunction for each 
^2 sector. To help ensure convergence, we then worked 
backwards, using the Sz — 74 solution for the Sz — 73 
guess and re-optimizing, retaining whichever wavefunc- 
tion gave the lowest energy before moving down to the 
next Sz sector. This sweeping procedure was especially 
helpful for resolving the minimum energies for Sz < 30. 



Many of the comparisons and insights we present in 
this section stem from the total energy results of our 
CPS ansatz, which are displayed in Table |l] and Fig- 
ure |3] As we described above, working with the CPS 
wavefunction in the variational Monte Carlo framework, 
it is simple to constrain the value of the total system's 
Sz spin, and so we are able to probe the lowest energy 
state in each Sz sector. As seen in Figure [3] the mini- 
mum energy as a function of Sz is nearly parabolic, as 
found in previous DMRG calculations [15^, and thus the 
CPS energies provide another wavefunction-based veri- 
fication of the qualitative correctness of the rotational 
band model (RBM). The agreement is not quantitative, 
however, and Figure |4] shows the deviation of the raw 
CPS energies when we try to fit them to the RBM form 
in Eq. ([3|. We see that except for the region near 1/3 
of the maximum magnetization, the differences between 
the CPS and RBM energies can be fit closely by a cubic 
correction, which is not surprising as cubic terms are the 
leading order terms neglected by the RBM. The sharp 
change in the deviations near 1/3 of saturation is more 
interesting, however, as it is responsible for creating the 
-Bsat/3 anomalies that cannot be predicted by the RBM. 
We will discuss these anomalies and the origins of the 



energy deviations responsible in Sections IV B and IV D 



TABLE I: Ground state energies of GPS and the RBM, in units of J, for the S — 5/2 Heisenberg model 
on the icosidodecahedron for different total Sz sectors of Hilbert space. The raw energies of the GPS 
wavefunction are given, as well as the energies produced by fits to the RBM using GPS energies, DMRG 
energies [IQ, and experimental magnetizations [14' . See Section IV A 
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-216.25 


-216.25 


-210.55 


-216.65 


38 


-61.23 


-62.52 


-58.15 


-62.77 


1 


-216.14 


-216.04 


-210.35 


-216.44 


39 


-52.91 


-54.43 


-50.13 


-54.67 


2 


-215.80 


-215.63 


-209.93 


-216.03 


40 


-44.39 


-46.14 


-41.90 


-46.36 


3 


-215.23 


-215.00 


-209.32 


-215.40 


41 


-35.68 


-37.63 


-33.47 


-37.85 


4 


-214.45 


-214.18 


-208.49 


-214.57 


42 


-26.78 


-28.92 


-24.83 


-29.13 


5 


-213.43 


-213.14 


-207.47 


-213.53 


43 


-17.67 


-20.00 


-15.99 


-20.20 


6 


-212.20 


-211.89 


-206.23 


-212.29 


44 


-8.37 


-10.87 


-6.94 


-11.06 


7 


-210.74 


-210.44 


-204.79 


-210.83 


45 


1.12 


-1.53 


2.31 


-1.71 


8 


-209.06 


-208.78 


-203.15 


-209.17 


46 


10.81 


8.01 


11.77 


7.84 


9 


-207.16 


-206.91 


-201.30 


-207.30 


47 


20.70 


17.76 


21.44 


17.60 


10 


-205.05 
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-199.24 
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62.16 
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-194.39 
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-194.84 


52 
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72.86 


69.52 


15 


-191.17 


-191.35 


-185.87 


-191.73 
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-187.76 
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-182.58 
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17 


-184.11 
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138.71 


141.34 
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59 
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150.95 


153.48 


150.92 
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69 
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284.76 


286.13 
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32 
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-106.71 
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-107.00 


70 
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35 
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FIG. 4: Deviation of the GPS energies from the fit to the 
RBM. Gircles represent the raw deviations, while lines repre- 
sent cubic fits in the ranges Sz £ [0, 19] and Sz G [30, 75]. See 
Section HVSI 



In the inset to Figure [3] we see that our CPS cal- 
culations produce superior variational energies as com- 
pared to DMRG. In addition to producing superior vari- 
ational energies, a fit of the CPS energies to the low- 
est band of the RBM produces band parameters {D = 
6.22,7 = 1.07) that resemble much more closely the 
band parameters fitted to experimental magnetization 
data {D ~ 6.23,7 = 1.07) [H] than those produced by a 
fit to the DMRG energies (D = 6.17,7 = 1.05) p] . 

In addition to RBM comparisons, we may compare the 
CPS singlet-triplet gap with the singlet-triplet gap mea- 
sured by neutron scattering. Using the value of J w 0.134 
meV [131 m], we find that our CPS calculations predict 
a gap of 0.015 meV, which is smaller than the RBM re- 
sult of 0.027 meV and significantly smaller than the 0.091 
meV gap derived from neutron scattering data |10| . Since 
both the Sz — and Sz — ^ energies are upper bounds, 
this suggests that the CPS Sz = Q ground-state energy, 
although an improvement over the DMRG energy, must 
still be relatively too high. This motivates further im- 
provements in the CPS ansatz, possibilities for which we 
mention in the conclusion. 



B. Properties 

Using the raw CPS energies or their fit to the rota- 
tional band model, we can evaluate the magnetization 
M, differential magnetization dM/dB, and heat capac- 
ity Cp of the S = 5/2 icosidodecahedron, as functions 
of the applied field B. Note that in our calculations of 
dM/dB and Cp we do not include the effects of excited 
states other than the lowest state in each spin-sector. In- 
deed, we show that including only the lowest spin-state 
contributions already produces much of the anomalous 
behavior common to spin systems built of corner shar- 
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FIG. 5: Total z magnetizations at different applied field 
strengths of the S = 5/2 icosidodecahedron. In the top panel 
we plot the low temperature limit of the magnetization curve 
derived from the GPS wavefunction energies. In the lower 
panel we plot finite temperature curves for the rotational band 
model using an ab initio parameterization based on the fit to 
the GPS energies (see Figure^, as well as the experimental 
results of Schnack et al [14]. See Section IV B 



ing triangles. Furthermore, in the case of dM/dB, the 
neglect of low-lying singlets is probably a reasonable ap- 
proximation at low temperatures, as such states do not 
contribute directly. 

To make quantitative predictions, we have taken |131 
[T3] the interaction strength as J/ks = 1.566iir and the 
spectroscopic splitting factor as g — 1.974. We begin 
by considering the magnetization curve at finite temper- 
atures, for which experimental results can be matched 
closely by the RBM [T3]. As may be expected by the 
similarity between the experimental and ab initio CPS 
fittings of the band model parameters discussed previ- 
ously, the magnetization curve obtained from the CPS- 
parameterized RBM also matches the experimental mag- 
netization curve closely, as seen in the bottom panel of 
Figure Is] More interesting, however, is the zero tem- 
perature limit of the magnetization curve, shown in the 
upper panel of Figure [51 where we see that the icosido- 
decahedron has anomalies in its magnetization staircase 



at field strengths near 1/3 of the saturation field strength 
Bsat = 17.7T. The staircase anomaly directly reflects the 
deviations in the ground-state energy levels, for each S"^, 
from the parabolic trend assumed by the RBM. We leave 
the discussion of the physical interpretation of the energy 
level deviations and corresponding staircase anomaly to 
section IfVDI Here we will show that this feature of the 
ground-state spectrum is sufficient to reproduce most of 
the unusual properties of the magnet near Bgat/S, with- 
out the need to explicitly consider other excited states. 

The differential susceptibility derived from the CPS 
energies is shown in Figure |6] We see that there is a 
sharp rise followed by a depression in the differential 
susceptibility, which in the case of the 0.5K results can 
be clearly associated with the staircase anomalies, which 
show up as gaps in the delta function progression of the 
OK dM/dB curve. Note that the area in the trough is 
greater than that in the peak, which, in conjunction with 
inhomogencities in the interactions [18] of the Mo72Fe3o 
compound that could smear the features together, may 
explain why only a broad trough is seen in experimental 
measurements [T3|TO]. In contrast to the CPS results, 
the RBM predicts only a very small dip in the dM/dB 
curve near Bsatl?>- 

As for the case of differential susceptibility, the heat 
capacity also shows a distinct feature near Bgat/^^ even 
when the low-lying excited states are ignored as in our 
CPS calculations. (Note that to the best of our knowl- 
edge, detailed measurements of the heat capacity are not 
yet available). As shown in Figure 7J the heat capacity 
derived from the CPS energies oscillates near the stair- 
case anomaly, whereas the heat capacity derived from 
the RBM shows only a small dip in this region. Note 
that the oscillations are present and essentially the same 
both when the CPS ground states are assumed to be non- 
degenerate and when they are assumed to have the same 
degeneracies as the corresponding states in the RBM. 
This offers reason to expect that the feature would be 
robust to the inclusion of additional excited states. 



TABLE IL Averages of the dot product Si- Sj/{\Si\\Sj\) for 
different choices of the sublattices (A,B,C) for sites i and j. 
The CPS results are for the wavefunction with zero total Sz, 
while the numbers for the classical Heisenberg model corre- 
spond to zero applied field. The abbreviation n.n. stands for 
nearest neighbor. See Section [IV C[ 
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FIG. 6: Differential susceptibility dM/dB as a function of the 
applied field strength in Teslas. Results for the CPS wave- 
function are shown both for the case when the ground state 
in each Sz sector is assumed to be non-degenerate (CPS) and 
when each Sz ground state is assumed to have the same de- 
generacy as the corresponding state in the RBM (CPS' ). For 
comparison, we also show the susceptibility derived from the 
first rotational band of the RBM with its experimentally de- 
rived parameterization [14]. Note that in the zero tempera- 
ture case the delta functions that make up the dM/dB curve 
have been scaled arbitrarily to show the number of magneti- 
zation levels ascended at each "step" , so for example the line 
just below 5T represents a magnetization change of 3, while 
most lines represent a change of 1. See Section [IV B[ 
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FIG. 7: Heat capacity as a function of the applied field 
strength in Teslas at T=0.4K. Results for the CPS wavefunc- 
tion are shown both for the case when the ground state in each 
Sz sector is assumed to be non-degenerate (CPS) and when 
each Sz ground state is assumed to have the same degeneracy 
as the corresponding state in the RBM (CPS^). For compar- 
ison, we also show the heat capacity derived from the first 
rotational band of the RBM with its experimentally derived 
parameterization [14]. See Section IV B 
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FIG. 8: The volumes of the parallelepiped defined by the spin 
triad's three vectors, given by the average of the scalar triple 
products S* ■ (S-* X S'') over all triangles with the sites ordered 
by sublattice {i £ A,j £ B,k £ C). Note that the individual 
triangles' triple products all had the same sign and that they 
deviated very little from the average. See Section [IV D| 
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FIG. 9: The expectation values of the product SlSiS^ of 
the three z components of the spin triad's vectors, averaged 
over all triangles {i,j,k}. Note that the individual triangles' 
products deviated very little from the average. See Section 

EYDl 



C. Spin correlations 



With the ground-state CPS wavefunction it is also 
possible to conipute the spin-spin correlation functions. 
These are shown in Table [ll] for the case of no external 
field. As the spin on the magnetic sites increases from 
1/2 to oo, the resulting ground-state is expected to be- 
come increasingly classical. The classical ground-state 
for corner-sharing triangles is well-known. Recall that 
the lattice is tripartite. Then all spins on sub-lattice A 
(and similarly for B and C) point in the same direction 
in the classical ground-state. The relative angle between 
the spins on sub-lattices A, B, and C is 120 degrees as is 
found in the classical ground-state for the Heisenberg tri- 
angle. Note that there are an infinity of classical ground- 
states, as the plane of the spins for sub-lattices A, B, and 
C can be rotated continuously. 

From our calculated correlation functions, the 
strongest correlations are naturally within the triangles. 
In the classical case, the spins are perfectly rotated from 
each other by 120 degrees, producing a dot product 
Si ■ Sj/{\Si\\Sj\) of -0.5 between nearest neighbor spins. 
Our CPS ansatz predicts that quantum fluctuations en- 
hance the expectation value of this dot product to -0.57. 
In doing so, the parallelity of spins on the same sub- 
lattice is disrupted. In the classical case we expect the 
spins to be perfectly parallel on the same sub-lattice, but 
quantum fluctuations reduce the average same sub-lattice 
dot product from 1.0 (the classical value) to 0.79. These 
values are the same for each sub-lattice, showing that at 
least by this metric, our ansatz preserves the equivalence 
of the different sub-lattices. 



D. Remnants of the Bsat/3 phase transition 

We now seek to provide some qualitative understand- 
ing of our wavefunction at different total Sz values. In 
doing so we will show that our wavefunction undergoes 
a change similar to a phase transition as the applied 
field is increased. In classical corner sharing triangle lat- 
tices, phase transitions are known to occur near Bsat/3 
between phases in which the spins on the different sub- 
lattices take on "Y" or "V" shaped configurations [37j . 
We must stress, however, that our numerical results do 
not allow us to fully distinguish whether our wavefunc- 
tion's phase transition behavior is a true property of the 
S = 5/2 icosidodecahedron or an artifact of our approx- 
imate ansatz, a point we discuss in some detail below. 

To characterize our wavefunction, we will focus on the 
behavior of the spin triads that make up each of the 
icosidodecahedron's triangles, an approach similar to the 
characterization of phases in the two dimensional infinite 
triangular and Kagome lattices. There, the phases are 
coplanar and described as either "Y" (or umbrella) , "V" , 
uud, or uuu [37] . The "Y" and "V" phases are so named 
because the shapes of these letters correspond to how the 
three spins are arranged in the plane (in the "V" case two 
of the spins are coUinear). We will see that our wave- 
function undergoes a sharp change between two states 
similar to the classical "Y" and "V" states, although our 
"Y" state is not coplanar and the spins in our "V" state 
may not be completely collinear (see Figures IS] and [9] for 
cartoons). 

To probe the character of our wavefunction's spin tri- 
ads, we have computed two expectation values. First, 
we have computed the scalar triple product of the three 
spin vectors of each triangle S^ ■ {S^ x S'^), which gives the 




FIG. 10: Evolution of a triad of classical spin vectors depicting the qualitative changes in our wavefunction with increasing 
applied field. The positions of the three vectors' endpoints on the S^ = 35/4 sphere are given by a red line with circles, a green 
line with squares, and a blue line with diamonds. The heavy black lines represent the spin vectors at i? = 0, the dashed grey 
lines represent the spin vectors just before the transition, and the dot-dot-dashed pink lines represent the spin vectors ju st afte r 
the transition. The thin black circle represents the intersection of the S^ — 35/4 sphere with the xz plane. See Section 



IV D 



volume of the parallelepiped that they define. For copla- 
nar or coUinear spins, the triple product will be zero, 
which should help us differentiate between these config- 
urations and others, such as a partially folded umbrella 
arrangement. In Figure Isj we plot the averages (over 
the twenty triangles) of the parallelepiped volumes for 
the ground state wavefunctions at different applied fields. 
We see that before Bsat/'i, the volume increases with field 
strength, which suggests the state may be a "folding um- 
brella" in which the initially 120° rotated spins gradu- 
ally close towards the z axis. However, near Bgat/^i the 
volume drops abruptly to zero and remains there for all 
higher field strengths. It is tempting to interpret this 
rapid drop as the remnants of what in classical 2D lat- 
tices would be a phase transition between non-coplanar 
and coplanar phases, although a recent study of the clas- 



sical triangular and Kagome lattices observes only copla- 
nar phases at low temperatures [37) . It appears that ei- 
ther quantum effects or errors inherent to our ansatz are 
stabilizing a non-coplanar arrangement for applied fields 
below Bsat/i- 

To further elucidate the qualitative nature of the states 
before and after Bgat/^^ we have also computed for each 
triangle the expectation values of the product of the three 
spins' 5*2 operators, SIS^S^. The average of these quan- 
tities over all triangles is shown in Figure |9] for different 
applied field strengths. We see that before Bsat/'i, the 
SzSzSz expectation values are negative, indicating that 
if the state is indeed a folding umbrella in this regime 
that one of the three spin vectors lies below the xy plane 
(this is the only way for three vectors with a non-negative 
net Sz to give a negative SzSzSz product). Thus it ap- 
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pears the state may be a tilted umbrella, in which the 
axis (or "handle" ) of the umbrella has been rotated away 
from vertical in such a way as to place one of the spokes 
below the xy plane. As with the triple product, our Sz 
product shows an abrupt change near Bsat/5, dropping 
rapidly to a value near the maximum-magnitude nega- 
tive of (—5/2)'^ characteristic of the uud configuration. 
It has been shown [I7J that the uud state makes a major 
contribution to the properties of the classical spin icosido- 
decahedron near Bsat/i- The analogous quantum states 
play a similar role in the 5 = 1/2 icosidodecahedron [T^ . 
and the same now appears to be true for S = 5/2. Upon 
increasing the field further, the Sz product rises smoothly 
to its maximum value at saturation, indicating that the 
three coplanar spin vectors are smoothly converting from 
an uud type configuration into the uuu configuration. 

As a final means to give a qualitative feel for our wave- 
function, we have constructed the classical spin triad that 
most closely matches the above expectation values. To 
do so, we have required that the triple products, Sz prod- 
ucts, and total Sz magnetizations match those given by 
our quantum wavefunction. In addition, we have arbi- 
trarily restricted one of the three vectors to the xz plane 
(our quantum expectation values are all rotationally in- 
variant about the z axis). Finally, in order to create a 
unique classical state, we have also required that the x 
and y components of the classical vectors each add to 
zero, as is the case for the classical ground state on a 
spin triangle [37] . These requirements give a unique evo- 
lution of the classical spin vectors with increasing applied 
field strength, which is depicted in Figure 10 Note that 



as the state approaches saturation, the quantum expec- 
tation values become incompatible with a classical spin 
state, and so we have only plotted the spin evolution up 
to the point at which compatibility fails. 

Unfortunately, we cannot rule out the possibility that 
the phase change behavior we observe may be an artifact 
of our approximate wavefunction. Our primary concern 
is that our initial guess is biased towards a particular col- 
oring of the icosidodecahedron lattice (see Section III C 
and Figure l2|, but this is not the only way to color the 
lattice and thus the initial guess does not possess all the 
correct symmetries. While it is possible that the opti- 
mization repairs this deficiency, it would be preferable 
to work with a wavefunction without this handicap. In 
future work it may be possible to use as the ansatz a lin- 
ear combination of CPS states with an initial guess taken 
such that each CPS in the combination is biased towards 
a different lattice coloring, thus removing fears of a col- 
oring bias in the overall ansatz. While our computer im- 
plementation is not currently capable of optimizing such 
an ansatz, we do not foresee any fundamental barriers to 
executing such an optimization in the future. 



V. CONCLUSIONS 

In this work we demonstrated that the correlator prod- 
uct state, a very simple ansatz designed for the treatment 
of strongly correlated spins, can successfully be used to 



model the quantum states of complex molecular mag- 
nets such as the Feao Keplerate system. The size of this 
system lies far outside the range of exact diagonalization. 
Our calculated variational energies are significantly lower 
than those previously obtained with the density matrix 
renormalization group and produce a fit to the rotational 
band model that is almost identical to that derived from 
experimental magnetization data. Furthermore, unlike 
the rotational band model, our ansatz is capable of pre- 
dicting anomalies in the differential susceptibility and 
heat capacity that are observed in frustrated magnetic 
systems near 1/3 of the saturation field. We have also 
analyzed a number of correlation functions among the 
spins, showing how the quantum state deviates from clas- 
sical behavior. Finally, we have shown how as a function 
of magnetic field, the quantum state appears to undergo 
a change reminiscent of phase transitions seen in classical 
2D corner sharing triangular lattices. 



In future research, more work is needed to clarify 
a number of aspects of this study. While our varia- 
tional energies are superior to previous theoretical treat- 
ments, neutron scattering data suggests there is still am- 
ple room for improvement, especially for small applied 
fields, as our predicted singlet-triplet gap is too small. 
In addition, it is not clear that the optimization of our 
ansatz fully preserves all symmetries of the icosidodeca- 
hedron, which makes definitive conclusions regarding the 
observed phase transition behavior difficult. To address 
these shortcomings, we have suggested that an ansatz 
consisting of a specially crafted linear combination of cor- 
relator product states be employed. The other obvious 
omission is a treatment of excited states, given that the 
presence of low-lying excitations is a key feature of frus- 
trated spin systems. Generalizing the methodology to 
model excited states within the Monte Carlo framework 
is not as straightforward as the generalization to linear 
combinations, but the critical importance of low-lying 
excitations makes it a highly desirable goal. 



The CPS family of states can naturally be applied 
to other magnetic systems as well as to more general 
non-spin electronic systems [10]. In the latter case, it 
is more advantageous to combine the correlators with a 
fermionic reference function. Correlators used in this way 
are formally the same as the Jastrow factors long studied 
in electronic structure. Jastrow factors are usually em- 
ployed to model "weak" correlations associated with the 
electron-electron cusp, while correlators have proven ef- 
fective at introducing strong correlations in the Hubbard 
model and some molecular systems. In the present study, 
we have shown that even the very complex correlations 
arising from magnetic frustration can be described effec- 
tively using correlators. Taken together, these findings 
motivate the use of correlator product states as a means 
to describe both weak and strong electron correlations 
simultaneously, a prospect that is under active investiga- 
tion. 
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